{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "baba5557",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from scipy.optimize import fsolve\n",
    "from scipy.stats import norm\n",
    "from matplotlib import pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8e55631a",
   "metadata": {},
   "source": [
    "# Identification of the parameters $\\{p, \\phi, \\alpha_0, \\alpha_1\\}$\n",
    "\n",
    "The model is\n",
    "$$\n",
    "\\max_{m_i} \\frac{(y_i-pm_i)^{1-\\sigma}}{1-\\sigma} + \\phi (1-\\exp(-\\alpha_1 - \\alpha_0 m_i))\n",
    "$$\n",
    "The FOC is \n",
    "$$\n",
    "p (y_i-pm_i)^{-\\sigma} = \\phi\\alpha_0 \\exp(-\\alpha_1 - \\alpha_0 m_i)\n",
    "$$\n",
    "This gives a solution $m_i = {\\cal M}(y_i,p|\\sigma, \\phi, \\alpha_0, \\alpha1)$\n",
    "\n",
    "\n",
    "### US\n",
    "For the US, we have $p=1$ and $\\sum_i \\omega_i y_i \\equiv y = 1$. Therefore, the moments restrictions  \n",
    "$$\n",
    "\\Psi_1 = \\frac{p \\sum_i \\omega_i m_i}{\\sum_i \\omega_i y_i}\n",
    "$$\n",
    "$$\n",
    "\\Psi_2 = \\sum_i (1-\\exp(-\\alpha_1 - \\alpha_0 m_i))\n",
    "$$\n",
    "$$\n",
    "\\Psi_3 = \\frac{(1-\\exp(-\\alpha_1 - \\alpha_0 m_{\\max}))}{(1-\\exp(-\\alpha_1 - \\alpha_0 m_{\\min}))}\n",
    "$$\n",
    "allows to identify $\\{\\phi,\\alpha_0,\\alpha_1\\}$\n",
    "\n",
    "\n",
    "### Europe\n",
    "\n",
    "For Europe, we keep $\\{\\phi,\\alpha_0,\\alpha_1\\}$ found in the previous step, and the moments restrictions\n",
    "$$\n",
    "\\Psi_1 = \\frac{p \\sum_i \\omega_i m_i}{\\sum_i \\omega_i y_i}\n",
    "$$\n",
    "$$\n",
    "\\Psi_2 = \\sum_i (1-\\exp(-\\alpha_1 - \\alpha_0 m_i))\n",
    "$$\n",
    "allows to identify $\\{p,\\alpha_1\\}$, given that a calibration of $\\sum_i \\omega_i y_i \\equiv y$\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1e3d32d8",
   "metadata": {},
   "source": [
    "## Functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "5f76a0af",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Function for finding the decision rule for m\n",
    "\n",
    "def medexp(m, *param_m):\n",
    "    p, y, sigma, alpha0, phi, alpha1 = param_m \n",
    "    return p*( (y - p*max(0,min(m,y/p)))**(-sigma) ) - alpha0*phi*np.exp(-( alpha1 + alpha0*max(0,min(m,y/p)) )) \n",
    "\n",
    "# Function for finding parameters (alpha0, phi, alpha1) using US moments' restrictions\n",
    "def hetero_calib(pd_c, *param_c):\n",
    "    \n",
    "    Psi1, Psi2, Psi3, p, yy, sigma = param_c \n",
    "    \n",
    "    msolv = np.zeros(yy.size)\n",
    "    \n",
    "    for ii in range(yy.size):\n",
    "        param_cm  = (p, yy[ii], sigma, pd_c[0], pd_c[1], pd_c[2])\n",
    "        m0        = .1\n",
    "        msolv[ii] = max(0,fsolve(medexp, m0, args=param_cm))\n",
    "    \n",
    "    return [Psi1 - np.mean(msolv), \n",
    "            Psi2 - np.mean(1-np.exp(-pd_c[2]-pd_c[0]*msolv)),\n",
    "            Psi3 - (1-np.exp(-pd_c[2]-pd_c[0]*msolv[-1]))/(1-np.exp(-pd_c[2]-pd_c[0]*msolv[0])) ] \n",
    "\n",
    "# Function for finding (p, alpha1) using European moments' restrictions\n",
    "def hetero_calib_eu(pd_s, *param_s):\n",
    "    \n",
    "    Psi1, Psi2, yy, sigma, phi, alpha0 = param_s \n",
    "    \n",
    "    msolv = np.zeros(yy.size)\n",
    "    \n",
    "    for ii in range(yy.size):\n",
    "        param_sm  = (pd_s[0], yy[ii], sigma, alpha0, phi, pd_s[1])\n",
    "        m0        = .1\n",
    "        msolv[ii] = max(0,fsolve(medexp, m0, args=param_sm))\n",
    "    \n",
    "    return [Psi1 - (pd_s[0]*np.mean(msolv))/(np.mean(yy)), \n",
    "            Psi2 - np.mean(1-np.exp(-pd_s[1]-alpha0*msolv))] \n",
    "\n",
    "# Function for finding the price, given all other parameters (alpha0, phi, alpha1)\n",
    "def hetero_new_price(p_new, *param_s):\n",
    "    \n",
    "    Psi1, yy, sigma, phi, alpha1, alpha0 = param_s \n",
    "    \n",
    "    msolv = np.zeros(yy.size)\n",
    "    \n",
    "    for ii in range(yy.size):\n",
    "        param_sn  = (p_new, yy[ii], sigma, alpha0, phi, alpha1)\n",
    "        m0        = .05\n",
    "        msolv[ii] = max(0,fsolve(medexp, m0, args=param_sn))\n",
    "    \n",
    "    return Psi1 - (p_new*np.mean(msolv))/(np.mean(yy)) "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e8bf25d4",
   "metadata": {},
   "source": [
    "## Calibrated parameter"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "49f43316",
   "metadata": {},
   "outputs": [],
   "source": [
    "sigma  = 2.0\n",
    "Ncount = 2"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bf29370a",
   "metadata": {},
   "source": [
    "## Aggregate moments\n",
    "\n",
    "Health inequalities: we use the gradient for the 4th quartile relative to the first. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "4acc8877",
   "metadata": {},
   "outputs": [],
   "source": [
    "grad_c  = [1.061, 1.276]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0365aed2",
   "metadata": {},
   "source": [
    "GDP share of health expenditures"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "ac3df66c",
   "metadata": {},
   "outputs": [],
   "source": [
    "s_c  = [0.0955, 0.1504]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0b855728",
   "metadata": {},
   "source": [
    "Income"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "20c3dd9c",
   "metadata": {},
   "outputs": [],
   "source": [
    "y_c  = [0.78, 1]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "aa65a7cb",
   "metadata": {},
   "source": [
    "For health status, we use the average health from the ADL definition in Table 4. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "be08d6d7",
   "metadata": {},
   "outputs": [],
   "source": [
    "h_c  = [0.931, 0.892]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "a1eab219",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0.0955, 0.931 , 1.061 ],\n",
       "       [0.1504, 0.892 , 1.276 ]])"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "psi = np.zeros((Ncount,3))\n",
    "for ii in range(Ncount):\n",
    "    psi[ii,:] = [s_c[ii],h_c[ii],grad_c[ii]]\n",
    "\n",
    "psi"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d01473d0",
   "metadata": {},
   "source": [
    "## Income inequalities: the US\n",
    "\n",
    "Table 5 gives the variance of the income process. If log income has mean 0 and variance equal to $\\sigma^2_{\\epsilon}$, we can solve for the quartiles (we use mid-points of quartiles, 12.5, 37.5, 62.5 and 87.5)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "48a95d99",
   "metadata": {},
   "source": [
    "We take the mean for EU in Table 5. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "e4046ef5",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[0.19111728571428574, 0.486429]"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sig_eu  = np.mean([0.256367, 0.094643, 0.214538, 0.237439, 0.169834, 0.094643, 0.270357])\n",
    "sig_c   = [sig_eu, 0.486429]\n",
    "sig_c"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0e5a0143",
   "metadata": {},
   "source": [
    "We find the percentiles, transform in levels and normalize to 1. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "a9f5396e",
   "metadata": {},
   "outputs": [],
   "source": [
    "qs   = [0.125,0.375,0.625,0.875]\n",
    "yy_c = np.zeros((Ncount,4))\n",
    "\n",
    "for ii in range(Ncount):\n",
    "    yy    = [norm(0,np.sqrt(sig_c[ii])).ppf(q) for q in qs]\n",
    "    yy    = np.exp(yy)\n",
    "    yy_c[ii,:] = y_c[ii]*yy/np.mean(yy)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f9017749",
   "metadata": {},
   "source": [
    "Here are relative incomes. Much more inequality in US"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "806c601a",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(4.975918497985015, 2.7340817164336597)"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "yy_c[-1,-1]/yy_c[-1,0], yy_c[0,-1]/yy_c[0,0]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0f2a9f4f",
   "metadata": {},
   "source": [
    "## Solution for the US"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "95c2f74f",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/loulou/anaconda3/lib/python3.8/site-packages/scipy/optimize/_minpack_py.py:175: RuntimeWarning: The iteration is not making good progress, as measured by the \n",
      "  improvement from the last ten iterations.\n",
      "  warnings.warn(msg, RuntimeWarning)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "array([7.51739835, 3.00535761, 1.46651755])"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "param_c   = (psi[1,0], psi[1,1], psi[1,2], 1, yy_c[1], sigma)\n",
    "pd_c0     = [6 , 2 , 2 ] \n",
    "calibsol  = fsolve(hetero_calib, pd_c0, args=param_c) \n",
    "calibsol"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8945fb87",
   "metadata": {},
   "source": [
    "We check the solution: to find the health gradient with the model "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "214a0e9a",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([0.        , 0.081768  , 0.18348814, 0.33634386]), 1.2759999999998748)"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "msolv = np.zeros(yy.size)\n",
    "\n",
    "for ii in range(yy.size):\n",
    "    param_m   = (1, yy_c[1,ii], sigma, calibsol[0], calibsol[1], calibsol[2])\n",
    "    m0        = .1\n",
    "    msolv[ii] = max(0,fsolve(medexp, m0, args=param_m))\n",
    "\n",
    "msolv, ( 1-np.exp(-(calibsol[2]+calibsol[0]*msolv[-1])) )/( 1-np.exp(-(calibsol[2]+calibsol[0]*msolv[0])) )"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1543e6dd",
   "metadata": {},
   "source": [
    "## Solution for the Europe"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "986c98f9",
   "metadata": {},
   "outputs": [],
   "source": [
    "#calibsol_eu = np.zeros((1,2))\n",
    "param_s           = (psi[0,0], psi[0,1], yy_c[0], sigma, calibsol[1], calibsol[0])\n",
    "pd_c0             = [0.9 , calibsol[2]]\n",
    "calibsol_eu       = fsolve(hetero_calib_eu, pd_c0, args=param_s)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "50d82726",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.56730391, 1.89073796])"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "calibsol_eu"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5d3a89f8",
   "metadata": {},
   "source": [
    "## Individual decisions and health gradient: Europe"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "dd23a4f1",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([0.01549786, 0.09418598, 0.16101862, 0.25451866]), 1.129474624558522)"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "msolv_eu = np.zeros(yy.size)\n",
    "for jj in range(yy.size):\n",
    "    param_m         = (calibsol_eu[0], yy_c[0,jj], sigma, calibsol[0], calibsol[1], calibsol_eu[1])\n",
    "    m0              = .1\n",
    "    msolv_eu[jj]    = max(0,fsolve(medexp, m0, args=param_m))\n",
    "grads_eu = ( 1-np.exp(-(calibsol_eu[1]+calibsol[0]*msolv_eu[-1])) )/(1-np.exp(-(calibsol_eu[1]+calibsol[0]*msolv_eu[0])))\n",
    "    \n",
    "msolv_eu, grads_eu\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "597eefdf",
   "metadata": {},
   "source": [
    "## Demand : price elasticity in the US"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "f3322eb8",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0.2654691129774273,\n",
       " 0.3251883988243632,\n",
       " 0.7497906679729276,\n",
       " 0.6677553518387016)"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "msolvp = np.zeros(yy.size)\n",
    "msolvm = np.zeros(yy.size)\n",
    "\n",
    "for ii in range(yy.size):\n",
    "    param_m   = (1.1, yy_c[1,ii], sigma, calibsol[0], calibsol[1], calibsol[2])\n",
    "    m0        = .1\n",
    "    msolvp[ii] = max(0,fsolve(medexp, m0, args=param_m))\n",
    "    param_m   = (0.9, yy_c[1,ii], sigma, calibsol[0], calibsol[1], calibsol[2])\n",
    "    m0        = .1\n",
    "    msolvm[ii] = max(0,fsolve(medexp, m0, args=param_m))\n",
    "\n",
    "dpop    = (1.1-1)/1\n",
    "share0  = s_c[1]\n",
    "share1p = 1.1*np.mean(msolvp)\n",
    "share1m = 0.9*np.mean(msolvm)\n",
    "level0  = s_c[1]\n",
    "level1p = np.mean(msolvp)\n",
    "level1m = np.mean(msolvm)\n",
    "\n",
    "((share1p-share0)/share0)/dpop, -((share1m-share0)/share0 )/dpop, ((level1m-level0)/level0)/dpop, -((level1p-level0)/level0 )/dpop\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "1ac7c9c0",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0.29532875590089525, 0.7087730099058146)"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "((((share1p-share0)/share0)/dpop)+(-((share1m-share0)/share0)/dpop))/2, ((((level1m-level0)/level0)/dpop)+(-((level1p-level0)/level0)/dpop))/2\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2db395ff",
   "metadata": {},
   "source": [
    "In the U.S., the price elasticity of the GDP share of health expenditures ($s_{US} = p_{US}m_{US}/y_{US}$) is equal to the price elasticity of the health expenditures ($p_{US}m_{US}$) becuase $y_{US}=1$.\n",
    "\n",
    "But, we can compute the price elasticity of the expenditures ($m_{US}$)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3c4bbc61",
   "metadata": {},
   "source": [
    "## Demand : income elasticity in the US"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "92110032",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0.15878235825745676,\n",
       " 0.29290743542285125,\n",
       " 1.1746605940832033,\n",
       " 1.2636166918805645)"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "msolvp = np.zeros(yy.size)\n",
    "msolvm = np.zeros(yy.size)\n",
    "\n",
    "for ii in range(yy.size):\n",
    "    param_m   = (1, 1.1*yy_c[1,ii], sigma, calibsol[0], calibsol[1], calibsol[2])\n",
    "    m0        = .1\n",
    "    msolvp[ii] = max(0,fsolve(medexp, m0, args=param_m))\n",
    "    param_m   = (1, 0.9*yy_c[1,ii], sigma, calibsol[0], calibsol[1], calibsol[2])\n",
    "    m0        = .1\n",
    "    msolvm[ii] = max(0,fsolve(medexp, m0, args=param_m))\n",
    "\n",
    "dyoy    = (1.1-1)/1\n",
    "share0  = s_c[1]\n",
    "share1p = np.mean(msolvp)/1.1\n",
    "share1m = np.mean(msolvm)/0.9\n",
    "level0  = s_c[1]\n",
    "level1p = np.mean(msolvp)\n",
    "level1m = np.mean(msolvm)\n",
    "\n",
    "((share1p-share0)/share0)/dyoy, -((share1m-share0)/share0)/dyoy, ((level1p-level0)/level0)/dyoy, -((level1m-level0)/level0)/dyoy "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "f964f2f9",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0.22584489684015402, 1.2191386429818838)"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "((((share1p-share0)/share0)/dyoy)+(-((share1m-share0)/share0)/dyoy))/2, ((((level1p-level0)/level0)/dyoy)+(-((level1m-level0)/level0)/dyoy))/2\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3af07c7a",
   "metadata": {},
   "source": [
    "## Demand : price elasticity in EU"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "ba949ad1",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(-0.06359513248237238,\n",
       " 0.014952998010533725,\n",
       " 1.0404498527973622,\n",
       " 0.9226845436459403)"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "msolv_eup = np.zeros(yy.size)\n",
    "msolv_eum = np.zeros(yy.size)\n",
    "\n",
    "for jj in range(yy.size):\n",
    "    param_m         = (calibsol_eu[0]*1.1, yy_c[0,jj], sigma, calibsol[0], calibsol[1], calibsol_eu[1])\n",
    "    m0              = .1\n",
    "    msolv_eup[jj]   = max(0,fsolve(medexp, m0, args=param_m))\n",
    "    param_m         = (calibsol_eu[0]*0.9, yy_c[0,jj], sigma, calibsol[0], calibsol[1], calibsol_eu[1])\n",
    "    m0              = .1\n",
    "    msolv_eum[jj]   = max(0,fsolve(medexp, m0, args=param_m))\n",
    "\n",
    "    \n",
    "share0  = s_c[0] \n",
    "share1p = calibsol_eu[0]*1.1*np.mean(msolv_eup)/np.mean(yy_c[0,:])\n",
    "share1m = calibsol_eu[0]*0.9*np.mean(msolv_eum)/np.mean(yy_c[0,:])\n",
    "\n",
    "level0  = s_c[0]/(calibsol_eu[0]/y_c[0]) \n",
    "level1p = np.mean(msolv_eup) \n",
    "level1m = np.mean(msolv_eum) \n",
    "\n",
    "((share1m-share0)/share0)/dpop, -((share1p-share0)/share0)/dpop, ((level1m-level0)/level0)/dpop, -((level1p-level0)/level0)/dpop \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "8e41e4a9",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(-0.024321067235919328, 0.9815671982216513)"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "((((share1m-share0)/share0)/dpop)+(-((share1p-share0)/share0)/dpop))/2, ((((level1m-level0)/level0)/dpop)+(-((level1p-level0)/level0)/dpop))/2\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bdaae432",
   "metadata": {},
   "source": [
    "## Demand : income elasticity in EU"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "ee03d8bb",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0.668744263187841, 0.8545471896810742, 1.7356186895066232, 1.769092470712967)"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "msolv_eup = np.zeros(yy.size)\n",
    "msolv_eum = np.zeros(yy.size)\n",
    "\n",
    "for jj in range(yy.size):\n",
    "    param_m         = (calibsol_eu[0], 1.1*yy_c[0,jj], sigma, calibsol[0], calibsol[1], calibsol_eu[1])\n",
    "    m0              = .1\n",
    "    msolv_eup[jj]   = max(0,fsolve(medexp, m0, args=param_m))\n",
    "    param_m         = (calibsol_eu[0], 0.9*yy_c[0,jj], sigma, calibsol[0], calibsol[1], calibsol_eu[1])\n",
    "    m0              = .1\n",
    "    msolv_eum[jj]   = max(0,fsolve(medexp, m0, args=param_m))\n",
    "\n",
    "share0  = s_c[0] \n",
    "share1p = calibsol_eu[0]*np.mean(msolv_eup)/(1.1*np.mean(yy_c[0,:]))\n",
    "share1m = calibsol_eu[0]*np.mean(msolv_eum)/(0.9*np.mean(yy_c[0,:]))\n",
    "    \n",
    "level0  = s_c[0]/(calibsol_eu[0]/y_c[0])    \n",
    "level1p = np.mean(msolv_eup) \n",
    "level1m = np.mean(msolv_eum)\n",
    "\n",
    "((share1p-share0)/share0)/dyoy, -((share1m-share0)/share0)/dyoy, ((level1p-level0)/level0)/dyoy, -((level1m-level0)/level0)/dyoy \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "1e6bf63d",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0.7616457264344576, 1.752355580109795)"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "((((share1p-share0)/share0)/dpop)+(-((share1m-share0)/share0)/dpop))/2, ((((level1p-level0)/level0)/dpop)+(-((level1m-level0)/level0)/dpop))/2\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cf13ddbd",
   "metadata": {},
   "source": [
    "# Forecasting of the GDP share of health expenditures\n",
    "\n",
    "Observation: the shift in GDP per capita\n",
    "\n",
    "Implication: the implicit shift in healthcare prices"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a5fa8267",
   "metadata": {},
   "source": [
    "### GDP per capita\n",
    "\n",
    "GDP share of health expenditures: OECD data\n",
    "\n",
    "Real GDP per capita: WB data (World Development Indicators)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "id": "2d32641c",
   "metadata": {},
   "outputs": [],
   "source": [
    "# 1970\n",
    "GDPpc_eu70 = [17887.10176179620, 26927.95903692180, 17519.42717022000, 15729.85644157820, 21354.25810368360, 24388.29821013090, 11431.59455874400]\n",
    "PMoY_eu70  = [5.713, 7.676666666666667,   5.184, float(\"NAN\"), 5.434246479956248,  5.506, 3.145]\n",
    "GDPpc70    = [np.nanmean(GDPpc_eu70), 25279.18879113590]\n",
    "PMoY70     = [np.nanmean(PMoY_eu70) , 6.232]\n",
    "# 1990\n",
    "GDPpc_eu90 = [29473.85438511120,   39295.24344249190,   28512.10870941340,   27479.57862963770,   31093.54486501490,   34570.30573399390,   18962.44605133250]\n",
    "PMoY_eu90  = [8.031, 8.038, 7.958, 7.015, 7.082, 7.256, 6.101]\n",
    "GDPpc90    = [np.nanmean(GDPpc_eu90) , 39278.59534978750]\n",
    "PMoY90     = [np.nanmean(PMoY_eu90) , 11.273]\n",
    "\n",
    "# 2000\n",
    "GDPpc_eu00 = [34476.20802831070, 49241.92800413450, 33583.85744208620, 32337.89674255510, 40440.67534426340, 41117.15243035039, 23928.67165633810]\n",
    "PMoY_eu00  = [9.827999999999999, 8.103999999999999, 9.541, 7.58, 7.056, 7.412, 6.8160]\n",
    "GDPpc00 = [np.nanmean(GDPpc_eu00) , 48689.03128985530]\n",
    "PMoY00  = [np.nanmean(PMoY_eu00) , 12.502]\n",
    "\n",
    "# 2017\n",
    "GDPpc_eu10 = [42622.40993268510, 55735.76490114850, 37678.92729753470, 31231.66444646460, 46978.44880120770, 52576.81028210510, 27213.45166906160]\n",
    "PMoY_eu10  = [11.272, 10.218, 11.458, 8.9010, 10.143, 10.916, 8.840999999999999]\n",
    "GDPpc10    = [np.nanmean(GDPpc_eu10) , 58387.77580835719]\n",
    "PMoY10     = [np.nanmean(PMoY_eu10) , 17.149999999999999]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "id": "4e128ceb",
   "metadata": {},
   "outputs": [],
   "source": [
    "psi0    = [i/100 for i in PMoY90]\n",
    "psi1    = [i/100 for i in PMoY10]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ba04b0a3",
   "metadata": {},
   "source": [
    "## Income heterogeneity by country (data)\n",
    "We rescale the income level and compute the quartiles in 1970 and 2017"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "id": "aab01b00",
   "metadata": {},
   "outputs": [],
   "source": [
    "y70_c = y_c*np.array(GDPpc70)/np.array(GDPpc00)\n",
    "y90_c = y_c*np.array(GDPpc90)/np.array(GDPpc00)\n",
    "y00_c = y_c*np.array(GDPpc00)/np.array(GDPpc00)\n",
    "y10_c = y_c*np.array(GDPpc10)/np.array(GDPpc00)\n",
    "\n",
    "yy70_c = np.zeros((Ncount,yy.size))\n",
    "yy90_c = np.zeros((Ncount,yy.size))\n",
    "yy00_c = np.zeros((Ncount,yy.size))\n",
    "yy10_c = np.zeros((Ncount,yy.size))\n",
    "\n",
    "for ii in range(Ncount):\n",
    "    yy    = [norm(0,np.sqrt(sig_c[ii])).ppf(q) for q in qs]\n",
    "    yy    = np.exp(yy)\n",
    "    yy70_c[ii,:] = y70_c[ii]*yy/np.mean(yy)\n",
    "    yy90_c[ii,:] = y90_c[ii]*yy/np.mean(yy)\n",
    "    yy00_c[ii,:] = y00_c[ii]*yy/np.mean(yy)\n",
    "    yy10_c[ii,:] = y10_c[ii]*yy/np.mean(yy)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "id": "de283a97",
   "metadata": {},
   "outputs": [],
   "source": [
    "prevsol0 = [0,0]\n",
    "prevsol1 = [0,0]\n",
    "\n",
    "#                  pm/y   income     sigma  phi          alpha_1         alpha_0\n",
    "param_s        = (psi0[0], yy90_c[0], sigma, calibsol[1], calibsol_eu[1], calibsol[0])\n",
    "pd_c0          = .1 #.3\n",
    "prevsol0[0]    = list(fsolve(hetero_new_price, pd_c0, args=param_s))\n",
    "param_s        = (psi0[1], yy90_c[1], sigma, calibsol[1], calibsol[2], calibsol[0])\n",
    "pd_c0          = .8\n",
    "prevsol0[1]    = list(fsolve(hetero_new_price, pd_c0, args=param_s))\n",
    "\n",
    "param_s        = (psi1[0], yy10_c[0], sigma, calibsol[1], calibsol_eu[1], calibsol[0])\n",
    "pd_c0          = .5#.5\n",
    "prevsol1[0]    = list(fsolve(hetero_new_price, pd_c0, args=param_s))\n",
    "param_s        = (psi1[1], yy10_c[1], sigma, calibsol[1], calibsol[2], calibsol[0])\n",
    "pd_c0          = 1.2\n",
    "prevsol1[1]    = list(fsolve(hetero_new_price, pd_c0, args=param_s))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "id": "43bea0f2",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "([[0.23604486640425065], [0.48068777532142404]],\n",
       " [[0.519900853169945], [1.3770087808513367]])"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "prevsol0, prevsol1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "id": "825acd69",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAD4CAYAAADlwTGnAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAATK0lEQVR4nO3dYYhdZ37f8e9vZYtdSFvZ9dgISVRmGULVQLRCaAV5s8TZImlLxwvd1iqxhTHIohZsIKVM0xd1+srdZuNiEBJyV0SmS4QhKR420wijrAkLlaPx1qu11hGeCnc968GeOKwT4xIj+98Xc9Tc3r3SPTOa1ezq+X7gcM95nv9z7/PA5f50js6dm6pCktSeT633BCRJ68MAkKRGGQCS1CgDQJIaZQBIUqPuWO8JrMQ999xT27dvX+9pSNLPlVdeeeUvqmpiuP3nKgC2b9/O3Nzcek9Dkn6uJPnfo9q9BCRJjTIAJKlRBoAkNcoAkKRG9QqAJPuSXE4yn2R6RH+SPNP1X0yya6h/Q5L/meRbA213J3kxyRvd4103vxxJUl9jAyDJBuAYsB/YARxMsmOobD8w2W2HgeND/V8FXh9qmwbOVdUkcK47liTdIn3OAPYA81V1pao+As4AU0M1U8Bztew8sCnJZoAkW4EvAf9lxJjT3f5p4MHVLUGStBp9AmAL8NbA8ULX1rfmPwP/BvhkaMx9VbUI0D3eO+rFkxxOMpdkbmlpqcd0JUl99AmAjGgb/hGBkTVJ/gnwblW9suKZXXuSqpNVtbuqdk9M/MQX2SRJq9Tnm8ALwLaB463A2z1r/hnwT5McAD4N/N0k/7Wqfh14J8nmqlrsLhe9u9pFSLeD7dN/tN5T0M+wN5/60po/Z58zgAvAZJL7k2wEHgJmhmpmgEe6u4H2Au9X1WJV/duq2lpV27txf9J9+F8bc6jbPwS8cLOLkST1N/YMoKquJjkKnAU2AKeq6lKSI13/CWAWOADMAx8Cj/Z47aeA55M8BvwQ+MrqliBJWo1efwyuqmZZ/pAfbDsxsF/AE2Oe4yXgpYHj94AH+k9VkrSW/CawJDXKAJCkRhkAktQoA0CSGmUASFKjDABJapQBIEmNMgAkqVEGgCQ1ygCQpEYZAJLUKANAkhplAEhSowwASWqUASBJjTIAJKlRBoAkNapXACTZl+Rykvkk0yP6k+SZrv9ikl1d+6eT/FmS7yW5lOS3B8Y8meRHSV7ttgNrtyxJ0jhjfxIyyQbgGPBFYAG4kGSmqn4wULYfmOy2zwPHu8e/AX61qj5IcifwnST/varOd+OerqrfWbvlSJL66nMGsAeYr6orVfURcAaYGqqZAp6rZeeBTUk2d8cfdDV3dlut1eQlSavXJwC2AG8NHC90bb1qkmxI8irwLvBiVb08UHe0u2R0Ksldo148yeEkc0nmlpaWekxXktRHnwDIiLbhf8Vft6aqPq6qncBWYE+SX+r6jwOfBXYCi8DXR714VZ2sqt1VtXtiYqLHdCVJffQJgAVg28DxVuDtldZU1Y+Bl4B93fE7XTh8AjzL8qUmSdIt0icALgCTSe5PshF4CJgZqpkBHunuBtoLvF9Vi0kmkmwCSPIZ4NeAP++ONw+M/zLw2s0tRZK0EmPvAqqqq0mOAmeBDcCpqrqU5EjXfwKYBQ4A88CHwKPd8M3A6e5Ook8Bz1fVt7q+ryXZyfKlojeBx9dqUZKk8cYGAEBVzbL8IT/YdmJgv4AnRoy7CHzuOs/58IpmKklaU34TWJIaZQBIUqMMAElqlAEgSY0yACSpUQaAJDXKAJCkRhkAktQoA0CSGmUASFKjDABJapQBIEmNMgAkqVEGgCQ1ygCQpEYZAJLUKANAkhrVKwCS7EtyOcl8kukR/UnyTNd/Mcmurv3TSf4syfeSXEry2wNj7k7yYpI3use71m5ZkqRxxgZA93u+x4D9wA7gYJIdQ2X7gcluOwwc79r/BvjVqvplYCewr/vReIBp4FxVTQLnumNJ0i3S5wxgDzBfVVeq6iPgDDA1VDMFPFfLzgObkmzujj/oau7sthoYc7rbPw08eBPrkCStUJ8A2AK8NXC80LX1qkmyIcmrwLvAi1X1cldzX1UtAnSP96549pKkVesTABnRVn1rqurjqtoJbAX2JPmllUwwyeEkc0nmlpaWVjJUknQDd/SoWQC2DRxvBd5eaU1V/TjJS8A+4DXgne4y0WKSzSyfIfyEqjoJnATYvXv3cPD0tn36j1Y7VLe5N5/60npPQVoXfc4ALgCTSe5PshF4CJgZqpkBHunuBtoLvN99sE8k2QSQ5DPArwF/PjDmULd/CHjh5pYiSVqJsWcAVXU1yVHgLLABOFVVl5Ic6fpPALPAAWAe+BB4tBu+GTjd3Un0KeD5qvpW1/cU8HySx4AfAl9Zu2VJksbpcwmIqppl+UN+sO3EwH4BT4wYdxH43HWe8z3ggZVMVpK0dvwmsCQ1ygCQpEYZAJLUKANAkhplAEhSowwASWqUASBJjTIAJKlRBoAkNcoAkKRGGQCS1CgDQJIaZQBIUqMMAElqlAEgSY0yACSpUQaAJDWqVwAk2ZfkcpL5JNMj+pPkma7/YpJdXfu2JN9O8nqSS0m+OjDmySQ/SvJqtx1Yu2VJksYZ+5OQ3e/5HgO+CCwAF5LMVNUPBsr2A5Pd9nngePd4FfjNqvpukr8DvJLkxYGxT1fV76zdciRJffU5A9gDzFfVlar6CDgDTA3VTAHP1bLzwKYkm6tqsaq+C1BVfw28DmxZw/lLklapTwBsAd4aOF7gJz/Ex9Yk2c7yD8S/PNB8tLtkdCrJXX0nLUm6eX0CICPaaiU1SX4B+APgN6rqr7rm48BngZ3AIvD1kS+eHE4yl2RuaWmpx3QlSX30CYAFYNvA8Vbg7b41Se5k+cP/m1X1h9cKquqdqvq4qj4BnmX5UtNPqKqTVbW7qnZPTEz0mK4kqY8+AXABmExyf5KNwEPAzFDNDPBIdzfQXuD9qlpMEuAbwOtV9buDA5JsHjj8MvDaqlchSVqxsXcBVdXVJEeBs8AG4FRVXUpypOs/AcwCB4B54EPg0W74rwAPA99P8mrX9ltVNQt8LclOli8VvQk8vkZrkiT1MDYAALoP7NmhthMD+wU8MWLcdxj9/wNU1cMrmqkkaU35TWBJapQBIEmNMgAkqVEGgCQ1ygCQpEYZAJLUKANAkhplAEhSowwASWqUASBJjTIAJKlRBoAkNcoAkKRGGQCS1CgDQJIaZQBIUqMMAElqVK8ASLIvyeUk80mmR/QnyTNd/8Uku7r2bUm+neT1JJeSfHVgzN1JXkzyRvd419otS5I0ztgASLIBOAbsB3YAB5PsGCrbD0x222HgeNd+FfjNqvqHwF7giYGx08C5qpoEznXHkqRbpM8ZwB5gvqquVNVHwBlgaqhmCniulp0HNiXZXFWLVfVdgKr6a+B1YMvAmNPd/mngwZtbiiRpJfoEwBbgrYHjBf72Q7x3TZLtwOeAl7um+6pqEaB7vLf3rCVJN61PAGREW62kJskvAH8A/EZV/VX/6UGSw0nmkswtLS2tZKgk6Qb6BMACsG3geCvwdt+aJHey/OH/zar6w4Gad5Js7mo2A++OevGqOllVu6tq98TERI/pSpL66BMAF4DJJPcn2Qg8BMwM1cwAj3R3A+0F3q+qxSQBvgG8XlW/O2LMoW7/EPDCqlchSVqxO8YVVNXVJEeBs8AG4FRVXUpypOs/AcwCB4B54EPg0W74rwAPA99P8mrX9ltVNQs8BTyf5DHgh8BX1mxVkqSxxgYAQPeBPTvUdmJgv4AnRoz7DqP/f4Cqeg94YCWTlSStHb8JLEmNMgAkqVEGgCQ1ygCQpEYZAJLUKANAkhplAEhSowwASWqUASBJjTIAJKlRBoAkNcoAkKRGGQCS1CgDQJIaZQBIUqMMAElqlAEgSY3qFQBJ9iW5nGQ+yfSI/iR5puu/mGTXQN+pJO8meW1ozJNJfpTk1W47cPPLkST1NTYAkmwAjgH7gR3AwSQ7hsr2A5Pddhg4PtD3e8C+6zz901W1s9tmr1MjSfop6HMGsAeYr6orVfURcAaYGqqZAp6rZeeBTUk2A1TVnwJ/uZaTliTdvD4BsAV4a+B4oWtbac0oR7tLRqeS3NWjXpK0RvoEQEa01Spqhh0HPgvsBBaBr4988eRwkrkkc0tLS2OeUpLUV58AWAC2DRxvBd5eRc3/p6reqaqPq+oT4FmWLzWNqjtZVburavfExESP6UqS+ugTABeAyST3J9kIPATMDNXMAI90dwPtBd6vqsUbPem1/yPofBl47Xq1kqS1d8e4gqq6muQocBbYAJyqqktJjnT9J4BZ4AAwD3wIPHptfJLfB74A3JNkAfj3VfUN4GtJdrJ8qehN4PG1W5YkaZyxAQDQ3aI5O9R2YmC/gCeuM/bgddof7j9NSdJa85vAktQoA0CSGmUASFKjDABJapQBIEmNMgAkqVEGgCQ1ygCQpEYZAJLUKANAkhplAEhSowwASWqUASBJjTIAJKlRBoAkNcoAkKRGGQCS1CgDQJIa1SsAkuxLcjnJfJLpEf1J8kzXfzHJroG+U0neTfLa0Ji7k7yY5I3u8a6bX44kqa+xAZBkA3AM2A/sAA4m2TFUth+Y7LbDwPGBvt8D9o146mngXFVNAue6Y0nSLdLnDGAPMF9VV6rqI+AMMDVUMwU8V8vOA5uSbAaoqj8F/nLE804Bp7v908CDq5i/JGmV+gTAFuCtgeOFrm2lNcPuq6pFgO7x3lFFSQ4nmUsyt7S01GO6kqQ++gRARrTVKmpWpapOVtXuqto9MTGxFk8pSaJfACwA2waOtwJvr6Jm2DvXLhN1j+/2mIskaY30CYALwGSS+5NsBB4CZoZqZoBHuruB9gLvX7u8cwMzwKFu/xDwwgrmLUm6SWMDoKquAkeBs8DrwPNVdSnJkSRHurJZ4AowDzwL/Ktr45P8PvA/gF9MspDksa7rKeCLSd4AvtgdS5JukTv6FFXVLMsf8oNtJwb2C3jiOmMPXqf9PeCB3jOVJK0pvwksSY0yACSpUQaAJDXKAJCkRhkAktQoA0CSGmUASFKjDABJapQBIEmNMgAkqVEGgCQ1ygCQpEYZAJLUKANAkhplAEhSowwASWqUASBJjeoVAEn2JbmcZD7J9Ij+JHmm67+YZNe4sUmeTPKjJK9224G1WZIkqY+xAZBkA3AM2A/sAA4m2TFUth+Y7LbDwPGeY5+uqp3dNosk6ZbpcwawB5ivqitV9RFwBpgaqpkCnqtl54FNSTb3HCtJWgd9AmAL8NbA8ULX1qdm3Nij3SWjU0nuGvXiSQ4nmUsyt7S01GO6kqQ++gRARrRVz5objT0OfBbYCSwCXx/14lV1sqp2V9XuiYmJHtOVJPVxR4+aBWDbwPFW4O2eNRuvN7aq3rnWmORZ4Fu9Zy1Juml9zgAuAJNJ7k+yEXgImBmqmQEe6e4G2gu8X1WLNxrb/R/BNV8GXrvJtUiSVmDsGUBVXU1yFDgLbABOVdWlJEe6/hPALHAAmAc+BB690djuqb+WZCfLl4TeBB5fw3VJksbocwmI7hbN2aG2EwP7BTzRd2zX/vCKZipJWlN+E1iSGmUASFKjDABJapQBIEmNMgAkqVEGgCQ1ygCQpEYZAJLUKANAkhplAEhSowwASWqUASBJjTIAJKlRBoAkNcoAkKRGGQCS1CgDQJIa1SsAkuxLcjnJfJLpEf1J8kzXfzHJrnFjk9yd5MUkb3SPd63NkiRJfYwNgCQbgGPAfmAHcDDJjqGy/cBktx0GjvcYOw2cq6pJ4Fx3LEm6RfqcAewB5qvqSlV9BJwBpoZqpoDnatl5YFOSzWPGTgGnu/3TwIM3txRJ0kr0+VH4LcBbA8cLwOd71GwZM/a+qloEqKrFJPeOevEkh1k+qwD4IMnlHnPWePcAf7Hek/hZkP+43jPQdfgeHXCT79N/MKqxTwBkRFv1rOkz9oaq6iRwciVjNF6Suaravd7zkK7H9+hPX59LQAvAtoHjrcDbPWtuNPad7jIR3eO7/actSbpZfQLgAjCZ5P4kG4GHgJmhmhngke5uoL3A+93lnRuNnQEOdfuHgBduci2SpBUYewmoqq4mOQqcBTYAp6rqUpIjXf8JYBY4AMwDHwKP3mhs99RPAc8neQz4IfCVNV2ZxvGymn7W+R79KUvVii7JS5JuE34TWJIaZQBIUqMMgNtcko+TvDqwTXftbya5Z6DuC0m+tX4zVauSbE/y2lDbk0n+dZK9SV7u3ruvJ3lynaZ5W+rzPQD9fPs/VbVzvSchrdJp4J9X1fe6Py3zi+s9oduJASDpZ9m9wLW/GPAx8IP1nc7txUtAt7/PDF0C+hfrPSFpBZ4GLif5b0keT/Lp9Z7Q7cQzgNvf9S4Bjbr/13uCtR6u976rqvoPSb4J/GPgXwIHgS/cqond7jwDaNd7wOBvMNyNf3hL62P4vQgD78eq+l9VdRx4APjlJH//Fs/vtmUAtOsl4GH4f7/b8OvAt9dzQmpTVX0ALCZ5AJZ/LArYB3wnyZeSXPujkpPAx8CP12WityG/CXybS/Ix8P2Bpj+uqukkf4/lH+75Ryz/1dY/Bqar6pN1mKYa1/1Q1DH+9kzgP1XVN5OcAXax/CdmrgL/rqrOrtM0bzsGgCQ1yktAktQoA0CSGmUASFKjDABJapQBIEmNMgAkqVEGgCQ16v8CGqsA7b6EwnIAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "columns = ('EU', 'US')\n",
    "Price0 = sum(prevsol0,[])\n",
    "Price1 = sum(prevsol1,[]) \n",
    "gPrice = [((Price1[ii]/Price0[ii])**(1/(2017-1990)))-1 for ii in np.arange(len(columns))]\n",
    "plt.bar(np.arange(len(columns)), gPrice)\n",
    "plt.xticks(np.arange(len(columns)), labels=columns)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "id": "9bd05d30",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[0.02967688121397627, 0.03974933924369717]"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "gPrice"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "23335f30",
   "metadata": {},
   "source": [
    "### Size of the error when price changes are omitted\n",
    "\n",
    "If prices change, then the model predicts exactly the GDP share of health expenditures (targeted data)\n",
    "\n",
    "$\\Rightarrow$ If prices do not change, then it exists a gap between model and data: this gap gives the error if price change are neglected "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "id": "a119d338",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "([8.11248145971843, 13.91943187305987],\n",
       " [7.3544285714285715, 11.273],\n",
       " [10.44032415091521, 15.510530852904951],\n",
       " [10.249857142857142, 17.15])"
      ]
     },
     "execution_count": 32,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "msolv0_90 = np.zeros((2,yy.size))\n",
    "pmoys0_90 = [0,0]\n",
    "msolv0_10 = np.zeros((2,yy.size))\n",
    "pmoys0_10 = [0,0]\n",
    "\n",
    "\n",
    "for jj in range(yy.size):\n",
    "    param_m        = (calibsol_eu[0], yy90_c[0,jj], sigma, calibsol[0], calibsol[1], calibsol_eu[1])\n",
    "    m0             = .1\n",
    "    msolv0_90[0,jj]= max(0,fsolve(medexp, m0, args=param_m))\n",
    "    param_m        = (calibsol_eu[0], yy10_c[0,jj], sigma, calibsol[0], calibsol[1], calibsol_eu[1])\n",
    "    m0             = .1\n",
    "    msolv0_10[0,jj]= max(0,fsolve(medexp, m0, args=param_m))\n",
    "\n",
    "for jj in range(yy.size):\n",
    "    param_m        = (1, yy90_c[1,jj], sigma, calibsol[0], calibsol[1], calibsol[2])\n",
    "    m0             = .1\n",
    "    msolv0_90[1,jj]= max(0,fsolve(medexp, m0, args=param_m))\n",
    "    param_m        = (1, yy10_c[1,jj], sigma, calibsol[0], calibsol[1], calibsol[2])\n",
    "    m0             = .1\n",
    "    msolv0_10[1,jj]= max(0,fsolve(medexp, m0, args=param_m))\n",
    "\n",
    "price_bench = [calibsol_eu[0],1] \n",
    "for ii in range(2):            \n",
    "    pmoys0_90[ii]= 100*price_bench[ii]*np.mean(msolv0_90[ii,:])/np.mean(yy90_c[ii,:])\n",
    "    pmoys0_10[ii]= 100*price_bench[ii]*np.mean(msolv0_10[ii,:])/np.mean(yy10_c[ii,:])\n",
    "\n",
    "pmoys0_90, PMoY90, pmoys0_10, PMoY10"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "83f75020",
   "metadata": {},
   "source": [
    "Without price change, the model over-estimate the GDP share of health expenditures in the 1990s"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "id": "2cdc6ddb",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[0.7580528882898578, 2.6464318730598695]"
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "[pmoys0_90[ii]-PMoY90[ii] for ii in np.arange(2)]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "267e1682",
   "metadata": {},
   "source": [
    "Without price changes, the model under-estimate the GDP share of health expenditures in the 2010s for the U.S."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "id": "ef97b436",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[0.19046700805806793, -1.6394691470950473]"
      ]
     },
     "execution_count": 34,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "[pmoys0_10[ii]-PMoY10[ii] for ii in np.arange(2)]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ef478f4c",
   "metadata": {},
   "source": [
    "# Forcasting long run change in GDP share of Health Expenditures"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7df75ef3",
   "metadata": {},
   "source": [
    "### Long run changes in GDP per capita and income dispersion"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "id": "cc424891",
   "metadata": {},
   "outputs": [],
   "source": [
    "gGDPpc    = [((GDPpc10[ii]/GDPpc90[ii])**(1/(2017-1990))-1) for ii in np.arange(2)]\n",
    "agGDPpc   = np.array(gGDPpc)\n",
    "agGDPpc10 = np.array(GDPpc10)\n",
    "\n",
    "horizon = 4\n",
    "fGDPpc  = np.zeros((2,horizon))\n",
    "fy_c    = np.zeros((2,horizon))\n",
    "fyy_c   = np.zeros((2,horizon,yy.size))\n",
    "\n",
    "for jj in range(2):\n",
    "    for ii in range(horizon):\n",
    "        fGDPpc[jj,ii] = (GDPpc10[jj]*(1+gGDPpc[jj])**(2030+(ii*10)-2017))/GDPpc00[jj]\n",
    "        fy_c[jj,ii] = y_c[jj]*fGDPpc[jj,ii]\n",
    "        yy    = [norm(0,np.sqrt(sig_c[jj])).ppf(q) for q in qs]\n",
    "        yy    = np.exp(yy)\n",
    "        fyy_c[jj,ii,:] = fy_c[jj,ii]*yy/np.mean(yy)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "id": "2d40c0d9",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[1.35719765, 1.53905893, 1.74528919, 1.97915382],\n",
       "       [1.45139469, 1.68093419, 1.94677557, 2.25466004]])"
      ]
     },
     "execution_count": 36,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "fGDPpc"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "id": "2b8441b3",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[0.01265430880932894, 0.014790792069827585]"
      ]
     },
     "execution_count": 37,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "gGDPpc"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a307dd87",
   "metadata": {},
   "source": [
    "### Long run changes in prices"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "id": "0f77bc8f",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0.76038411, 1.01869143, 1.36474739, 1.82836075],\n",
       "       [2.28564722, 3.3751706 , 4.9840485 , 7.35984708]])"
      ]
     },
     "execution_count": 38,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "fPrice  = np.zeros((2,horizon))\n",
    "\n",
    "for jj in range(2):\n",
    "    for ii in range(horizon):\n",
    "        fPrice[jj,ii] = (Price1[jj]*(1+gPrice[jj])**(2030+(ii*10)-2017))\n",
    "        \n",
    "fPrice"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7e29d255",
   "metadata": {},
   "source": [
    "Europe with changes in price"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "id": "bbfbdab0",
   "metadata": {},
   "outputs": [],
   "source": [
    "xx     = np.zeros((horizon,yy.size))\n",
    "yy_eu  = np.zeros(horizon)\n",
    "yy_eu0 = np.zeros(horizon)\n",
    "yy_us  = np.zeros(horizon)\n",
    "yy_us0 = np.zeros(horizon)\n",
    "\n",
    "for ii in range(horizon):\n",
    "    for jj in range(yy.size):\n",
    "        param_m          = (fPrice[0,ii], fyy_c[0,ii,jj], sigma, calibsol[0], calibsol[1], calibsol_eu[1])\n",
    "        m0               = .05#.01\n",
    "        xx[ii,jj] = max(0,fsolve(medexp, m0, args=param_m)) \n",
    "    \n",
    "    yy_eu[ii] = 100*fPrice[0,ii]*np.mean(xx[ii,:])/np.mean(fyy_c[0,ii,:])\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c852a7b4",
   "metadata": {},
   "source": [
    "Europe without price changes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "id": "ce3c0840",
   "metadata": {},
   "outputs": [],
   "source": [
    "for ii in range(horizon):\n",
    "    for jj in range(yy.size):\n",
    "        param_m          = (Price1[0], fyy_c[0,ii,jj], sigma, calibsol[0], calibsol[1], calibsol_eu[1])\n",
    "        m0               = .01#.005\n",
    "        xx[ii,jj] = max(0,fsolve(medexp, m0, args=param_m)) \n",
    "    \n",
    "    yy_eu0[ii] = 100*Price1[0]*np.mean(xx[ii,:])/np.mean(fyy_c[0,ii,:])\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "99c5497f",
   "metadata": {},
   "source": [
    "The U.S. with price changes "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "id": "ce7da493",
   "metadata": {},
   "outputs": [],
   "source": [
    "for ii in range(horizon):\n",
    "    for jj in range(yy.size):\n",
    "        param_m          = (fPrice[1,ii], fyy_c[1,ii,jj], sigma, calibsol[0], calibsol[1], calibsol[2])\n",
    "        m0               = .01#.005\n",
    "        xx[ii,jj] = max(0,fsolve(medexp, m0, args=param_m)) \n",
    "    \n",
    "    yy_us[ii] = 100*fPrice[1,ii]*np.mean(xx[ii,:])/np.mean(fyy_c[1,ii,:])\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a2079c7d",
   "metadata": {},
   "source": [
    "The U.S. without price changes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "id": "58fb5606",
   "metadata": {},
   "outputs": [],
   "source": [
    "xx = np.zeros((horizon,yy.size))\n",
    "for ii in range(horizon):\n",
    "    for jj in range(yy.size):\n",
    "        param_m          = (Price1[1], fyy_c[1,ii,jj], sigma, calibsol[0], calibsol[1], calibsol[2])\n",
    "        m0               = .01#.005\n",
    "        xx[ii,jj] = max(0,fsolve(medexp, m0, args=param_m)) \n",
    "    \n",
    "    yy_us0[ii] = 100*Price1[1]*np.mean(xx[ii,:])/np.mean(fyy_c[1,ii,:])\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "id": "8003a56b",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "The PostScript backend does not support transparency; partially transparent artists will be rendered opaque.\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAD4CAYAAAD2FnFTAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAA6yElEQVR4nO3deXhUVbbw/+/KxDwIQaYAgRAigwgyqBdEFIQINFFsEW3tdmgRFK+t9CBqi6/989dIe9u376VbRfQqraDcayCImIC0tgRHUFQEQhgCxIQhTCEkIdN6/6hKdSWpykRVKqmsz/PkqaqTc04tDpWzap+9z16iqhhjjDGVhQQ6AGOMMY2TJQhjjDEeWYIwxhjjkSUIY4wxHlmCMMYY41FYoAPwpcjISI2Ojg50GMYY02Rs27YtR1W7ePpdUCWI6Ohotm7dGugwjDGmyRCRg95+Z5eYjDHGeGQJwhhjjEeWIIwxxngUVH0QnhQXF5OZmUlhYWGgQwlKLVu2JCoqivDw8ECHYozxsaBPEJmZmbRr147o6GhEJNDhBBVV5cSJE2RmZtK3b99Ah2OM8TG/JQgR6QUsB7oBZcBSVf2LiLwDxDlX6wicVtVhHrbPAM4CpUCJqo6sTxyFhYWWHPxEROjcuTPHjx8PdCjGNDtjxowhJyenyvLIyEi2bNnik/fwZwuiBJivql+LSDtgm4hsVNVby1cQkf8AzlSzj2tVteoRqCNLDv5jx9aYwPCUHKpbXh9+SxCqmg1kO5+fFZFdQE9gJ4A4ziwzgev8FYMxxpj6a5BRTCISDQwHvnBbfDVwVFXTvWymwAYR2SYis/0cojHGNAnnzp1j/fr1DfJefu+kFpG2wLvAr1Q11+1XtwErq9l0jKpmicjFwEYR2a2qn3jY/2xgNkDv3r19GLkxxjQOeXl5fPzxxyQnJ/PJJ59w/vz5Bnlfv7YgRCQcR3J4S1UT3ZaHATOAd7xtq6pZzsdjwGpgtJf1lqrqSFUd2aWLx+lEGoXQ0FCGDRvm+lm0aBEZGRkMGTKkwnpPP/00zz//fI37W716NSLC7t27/RUyycnJxMXF0b9/fxYtWuS39zHGVJWXl8fatWt54IEHuOqqq5g/fz7ffvstM2fO5K233mqQGPw5ikmAV4FdqvrnSr+eCOxW1Uwv27YBQpx9F22AScAz/orVXWpqKqtWrSInJ4fIyEhmzpzJ2LFjL3i/rVq1Yvv27RWWZWRk1Ht/K1euZOTIkbz99ts8/fTTFxRbaWkpWVlZ9OrVq8KyBx98kI0bNxIVFcWoUaOYPn06gwYNuqD3MsZ4l5ubyz/+8Q+Sk5NJTU2luLiYrl27MmvWLOLj4xk+fDghIY7v9ZGRkV5HMfmKPy8xjQHuBL4Xke3OZY+r6npgFpUuL4lID2CZqk4BugKrnSNkwoAVqprsx1gBR3JYtmwZRUVFgGM0wLJlywB8kiR8JS8vj3/+859s3LiRW265pUKCmDVrFqpKRkYGR44c4W9/+xtTp06tdn/z588nJiaGhx56yLXsyy+/pH///vTr18+136SkJEsQxvjYmTNn2LRpEykpKWzZsoXi4mK6d+/Oz372M+Lj47nssstcScGdr4ayVsefo5hSAY9jIFX1Lg/LsoApzuf7gct8HdPy5cs5eNDrxIWkp6dTUlJSYVlRURFLly7lo48+8rhNnz59+PnPf17jexcUFDBs2DDX6wULFnDFFVfULvBK1qxZw8SJExk6dCht2rTh66+/5vLLLwfg22+/5cYbb+Sdd94hNTWVRx99tEqCmDhxIkeOHAEcN7vt2rWLQYMG0bt3bxISEgD48ccfK7QooqKi+OKLLzDGXLjTp0/z4YcfkpKSwmeffUZxcTE9e/bkzjvvJD4+nqFDhzaKIeRBfyd1XVRODjUtrwtPl5i8JauaPhgrV65k9mzHwK6ZM2eycuVKLr/8cgoKCsjJyWHhwoUADBo0iFOnTlXZ/sMPPwTg008/5be//S35+fm0bNmywjqqWue4jDHenTx5kk2bNpGcnMznn39OSUkJUVFR/PznPyc+Pp5LL7200f2NNasEUdM3/X//93/3ek3v97//vc/j6dy5c5UT+MmTJ6udtuLEiRN8+eWXJCY6+vxvvfVWrrnmGhYvXsyOHTuIjY11ney//vprLrusakOsvAVx+PBh2rdvz8iRjpvUn332WVcLIioqisOHD7u2yczMpEePHhf2DzammTl58iQbN24kOTmZL774gtLSUnr37s0999zD5MmTGTx4cKNLCu6aVYKoycyZMyv0QQBEREQwc+ZMv7xf27Zt6d69O5s2bWLChAmcPHmS5ORkHn74YQAmTJjA8uXL6dmzp2ub//3f/2XKlCm0aNECgL59+9KtWzdSU1NJS0vj0KFDFBYWUlpaysKFC1m8eHGV933vvfe4/vrrefXVV71e5ho1ahTp6ekcOHCAnj178vbbb7NixQo/HAVjgktOTo4rKXz55ZeUlZURHR3Nfffdx+TJkxk4cGCjTgruLEG4Ke+IXrVqFSdOnKBz584+G8VUuQ8iPj6eRYsWsXz5ch588EHmz58PwMKFC4mJiaGsrIy9e/fSqVOnCvtZuXIl3333He6lVU+cOMGKFSsICwvjZz/7GePHjyc3N5fHH3+cMWPGeIzn2WefrbYPJCwsjCVLljB58mRKS0u55557GDx4cP0PgDFB7Pjx42zYsIHk5GS2bt1KWVkZffv25f777yc+Pp64uLgmkxTciadrzU3VyJEjtXLJ0V27djFw4MAARVR/O3bs4LXXXuPPf648Qti7cePG8corrxAXF1fzyj7UVI+xMRfi6NGjbNiwgZSUFLZu3YqqEhMTQ3x8PPHx8cTGxjaJpCAi27xNhmotiEZqyJAhdUoOAPv27SM2NtZPERljjhw5QkpKCikpKXz99deoKrGxscybN4/JkycH3d+fJYgg8uOPPwY6BGOCTnZ2NikpKSQnJ/PNN98AMGDAAB566CHi4+OJiYkJcIT+YwnCGGMq+fHHH11J4dtvvwXgkksu4Ve/+hWTJ0923UAa7CxBGGMMcPjwYVdS+P777wEYPHgwjz76KJMnT64wMKS5sARhjGm2Dh06RHJyMsnJyfzwww+Ao/9v/vz5xMfHN/sZoi1BGGOalYyMDFdLYefOnQAMHTqU3/zmN0yePLnCFDPNnSUIY0zQ279/P8nJyaSkpLimyB82bBiPPfYYkyZNqnAzqvkXSxDGmKC0b98+PvjgA1JSUtizZw8Aw4cPZ8GCBUyePJnu3bsHOMLGzxKEMSZopKenu/oU9u7di4hw+eWX88QTTzBp0iS6desW6BCbFEsQxpgmS1XZs2eP6/LRvn37EBFGjhzJ73//e66//nq6du0a6DCbLL+WHDUONZUWffbZZxk8eDBDhw5l2LBhta674O+yo1Zy1DRGqsru3bt54YUXiI+PZ/r06bz00kt07tyZp556is2bN/Pmm29yxx13WHK4QNaCcDNmzBiv0337q3rTZ599xrp16/j6669p0aIFOTk5FWaTrY6vyo5ayVHTGFT395eamsrOnTtdo48OHjxISEgIV1xxBXfddRfXX3+9T0ttGgdLEG48fTirW+4L2dnZREZGuqbvru2H3FvZUSs5apqq6v7+rr/+eg4fPkxoaChXXnkl9957L9dff32V2Y6NbzWrBPHss8/W+3LMnXfe6XH5JZdcwhNPPFHvmCZNmsQzzzzDgAEDmDhxoqsAUE28lR21kqMmGPXp04f777+fCRMmWFJoQM0qQQSKtyl/RYS2bduybds2Nm/ezEcffcStt97KokWLuOuuu6rdp6eyowMHDrSSo6ZJysrKqvb3r776agNFYtz5LUGISC9gOdANKAOWqupfRORp4D7guHPVx1V1vYft44G/AKHAMlW94F7Smr7pV1dH4e9//3u937em0qKhoaGMHz+e8ePHc+mll/LGG29UmyC8lR2dOXOmlRw1TUZeXh4pKSkkJSXx5ZdfBjoc44E/RzGVAPNVdSBwJfCgiJRfwH5BVYc5fzwlh1Dgr8ANwCDgNrdtmxz30qKAq7To2LFjSUtLIz093bXu9u3b6dOnD+AoOeppCm9vZUeTk5NdJUfPnTvHwoULeeSRR6ps/95779GxY0c2bNjA4cOH2bFjBzt27HAlB6hYcrSoqIi3336b6dOn+/S4mOanpKSEf/7zn8yfP58xY8bw+OOPc+TIEebNmxfo0IwHfmtBqGo2kO18flZEdgG1vZ99NLBXVfcDiMjbQAKw0x+xlouMjPQ6iuJCeSstum3bNh566CFOnz5NWFgY/fv3Z+nSpV5LjoL3sqMjRoywkqOm0Snv51qzZg3vv/8+OTk5dOzYkZtuuomEhASGDRuGiLBy5Uq//f2Z+mmQkqMiEg18AgwBHgXuAnKBrThaGacqrf9TIF5Vf+l8fSdwhapW+ZohIrOB2QC9e/cecfDgwQq/b6rlMK3kqGnqjhw5wtq1a1m7di3p6emEh4czfvx4brzxRsaNG0dERESgQzQEuOSoiLQF3gV+paq5IvIi8AdAnY//AdxTeTMPu/KYyVR1KbAUHDWpfRV3oFnJUdMU5eXlsXHjRpKSkvj8889RVYYPH87TTz/NDTfcQMeOHQMdoqkDvyYIEQnHkRzeUtVEAFU96vb7V4B1HjbNBNzn3I0Cqh/mYKzkqAmIkpISPv30U5KSkvjwww8pLCykV69ePPjgg0yfPt3Vp2aaHn+OYhLgVWCXqv7ZbXl3Z/8EwE3ADg+bfwXEikhf4EdgFnC7v2I1xtTd7t27WbNmDevWreP48eO0b9+ehIQEEhISuPzyy21YdBDwZwtiDHAn8L2IbHcuexzHiKRhOC4ZZQD3A4hIDxzDWaeoaomIzANScAxzfU1Vf/BjrMaYWjh69CjvvfceSUlJ7Nmzh/DwcK655hoSEhIYP3689SsEGX+OYkrFc19ClWGtzvWzgClur9d7W9cY03DOnTvn6lf47LPPUFWGDRvGU089xZQpU7jooosCHaLxE7uT2hhTRWlpKZ999hlJSUls3LiRgoICoqKimDt3LgkJCRWGWJvgZQnCGOOye/dukpKSWLduHceOHaNdu3ZMnz6d6dOnM2LECOtXaGYsQRjTzB07dox169axZs0a0tLSCAsLY9y4cSQkJHDttde67tg3zY8lCGOaofz8fDZu3MjatWv59NNPKSsrY+jQofz+979nypQpNmOqASxBGNNslJaW8sUXX5CUlMSGDRvIz8+nZ8+e3H///UyfPt1V+8OYcpYgGkhoaCiXXnqp6/WsWbOYNWsW06ZNY8eOf90K8vTTT9O2bVt+/etfV7u/1atXM2PGDHbt2sUll1zil5iTk5N5+OGHKS0t5Ze//CWPPfaYX97H+NeePXtISkrivffe4+jRo7Rr146pU6eSkJDAiBEjCAmxysPGM0sQHpw/f56dO3cyaNAgn11/bdWqFdu3b6+wLCMjo97781W5UW+s5GjTdvz4cdatW0dSUhK7du0iLCyMq6++mgULFnDddddZv4KpFUsQHmRkZHDmzBkyMjIafOK72vBWbhTqV3LUEys52vQUFBTw4YcfkpSUxJYtWygrK2PIkCE8+eSTTJ061foVTJ01qwSRnp5OXl6e19+fOXOmwuvs7Gyysx2zgnTo0MHjNm3btq3VBHkFBQUMGzbM9XrBggXVTrddHW/lRoFalRy9+uqrOXv2bJX9Pv/880ycOBGwkqNNRVlZGV988QVr164lOTmZ/Px8evTowezZs5k+fToxMTGBDtE0Yc0qQdSkXbt2FBYWUlxc7FoWHh5epRxnfXi6xFR5avJyNY0191Ru9PLLL6egoKBWJUc3b95cY7xWcrRxS09Pd/UrHDlyhDZt2nDDDTeQkJDAqFGjrF/B+ESzShC1+aaflpZGdnY2ISEhlJWVERkZ6bfLTDWVIvXEW7nRxYsXs2PHjlqVHK1NC8JKjjY+OTk5rFu3jrVr1/LDDz8QGhrK2LFj+d3vfsd1113nky8yxrhrVgmiNoqLi+nRowc9evQgKyuLoqIiv72XeynSCRMmuEqRPvzww4Cj5Ojy5cvp2fNfhfi8lRtNTU0lLS3NVXK0tLSUhQsXsnjx4irvW5sWhHvJ0Z49e/L222+zYsUKH/3LTW0VFhayadMmkpKSSE1NpbS0lMGDB/P4448zdepUq7Zm/MoSRCVDhgxxPR8wYIDP9lu5DyI+Pp5FixZ5LUXqreSot3KjK1asICwsrNYlR2tiJUcDp6ysjK+++oqkpCSSk5M5d+4c3bt359577yUhIYH+/fsHOkTTTFiCaCClpaUelw8aNIiPPvqoyvKdO3dy880306pVqwrLP/74Y6/vUV5y9LnnnrugWMtNmTKFKVOm1Lyi8Yl9+/a5+hWysrJo3bo18fHxJCQkMHr0aOtXMA3OEkQjZSVHm4cTJ07w/vvvk5SUxI4dOwgJCWHs2LHMnz+fCRMmVPmCYExDsgQRRKzkaOMyZswYcnJyqizv3LkzTz75JElJSWzevJnS0lIGDRrEggULmDp1Kl26dAlAtMZUZQnCGD/xlBzA0Wp45JFH6Nq1K/fccw8JCQnW8jONkiUIYwLg9ddfZ/To0YSGhgY6FGO8sgRhjI+dO3eOlJSUate56qqrGigaY+rPEoQxPqCqbN26lcTERNeUF8Y0dX5LECLSC1gOdAPKgKWq+hcR+RPwE6AI2AfcraqnPWyfAZwFSoESVR3pr1iNqa/s7GzWrFlDYmIihw4donXr1kyZMoUZM2Zw++23Bzo8Yy6IP1sQJcB8Vf1aRNoB20RkI7ARWKCqJSLyHLAA+J2XfVyrqp57+owJkPK7mxMTE9myZQuqyujRo3nwwQeZNGkSrVu3BiAyMtJjR7Xd/WyaCr8lCFXNBrKdz8+KyC6gp6pucFvtc+Cn/orBGF9RVb7//nveffdd1q9fT25uLj169OCBBx7gpptuqjDzbbktW7YEIFJjfKdB+iBEJBoYDlSeL/oe4B0vmymwQUQUeFlVl3rZ92xgNkDv3r19Eq8x5XJyckhKSmL16tWkp6fTokULJk2axM0338wVV1xhdzeboOb3T7eItAXeBX6lqrluy5/AcRnqLS+bjlHVy4EbgAdFZJynlVR1qaqOVNWRvrrBKC0tjQULFpCWluaT/WVkZFSY4wkcpUWff/55AJ599lkGDx7M0KFDGTZsWK3rLqxevRoRYffu3T6Js7Lk5GTi4uLo378/ixYt8st7NEZFRUVs3LiROXPmMG7cOBYvXkybNm145pln2LJlC88//zxXXXWVJQcT9PzaghCRcBzJ4S1VTXRb/gtgGjBBPRUeAFQ1y/l4TERWA6OBT/wZLziSw5IlSygqKmLJkiXMmzfPr1XlPvvsM9atW8fXX39NixYtyMnJqfUMsv4sO9ocS47u3r2bxMRE3nvvPU6ePEmXLl24++67mTFjhhXeMc2SP0cxCfAqsEtV/+y2PB5Hp/Q1qupxLKCItAFCnH0XbYBJwDP+irWce3IAGiRJZGdnExkZ6Zq+u7YdmN7KjlrJ0bo5ffo069atIzExkR9++IHw8HCuu+46ZsyYwdixYwkLs5Hgpvny56d/DHAn8L2IbHcuexz4T6AFsNFZoexzVZ0jIj2AZao6BegKrHb+PgxYoarJFxrQO++8Q2Zmpsff5efn8+OPP1appFZUVMQLL7xAz549XaNT3EVFRXHrrbfWO6ZJkybxzDPPMGDAACZOnOgqAFQTb2VHreRozUpLS0lNTSUxMZFNmzZRXFzMwIEDeeKJJ5g2bZrVbjbGyZ+jmFIBTzUq13tZPwuY4ny+H6haCs2Pjhw54rHMpjMejhw54vo2XVfeSnWKCG3btmXbtm1s3ryZjz76iFtvvZVFixZx1113VbtPT2VHBw4caCVHq3HgwAESExNZs2YNx44do2PHjsyaNYubb76ZgQMHBjo8YxqdZtV+ru6bfuXLS+4iIiIu6DJTTaVFQ0NDGT9+POPHj+fSSy/ljTfeqDZBeCs7OnPmTCs5WkleXh4ffPAB7777Lt988w0hISGMGzeOJ598kmuvvZaIiIhAh2hMo9WsEkR14uLimDdvXpUkcaHJAaovLZqWlkZISIhrNs/t27fTp08fwHPJUfBedjQ5OdlKjuKoyFaeQDds2EBBQQH9+vXj17/+NQkJCVx88cWBDtGYJsEShJvKScIXyaGct9Ki27Zt46GHHuL06dOEhYXRv39/li5d6rXkKHgvOzpixIhmXXI0MzOTNWvWsHr1ajIzM2nbti3Tp09nxowZXHbZZU3+EpkxDU28XXdvikaOHKlbt26tsGzXrl11vr6clpbG66+/zl133eXXIa7V2bFjB6+99lqdqsqVlxxt6Jjrc4x9paCggA0bNpCYmMjnn3+OiHDVVVcxY8YMJk6caBXZjKmBiGzzNtedJYgg0rNnTw4fPtzgN3A19DFWVbZv305iYiLr168nLy+PqKgoZsyYwY033ljlkpwxxrvqEoRdYgoiwV5y9OjRo65pL/bv30+rVq2YPHkyM2bMYNSoUXZnszE+ZgnCNGpFRUX84x//IDExkc2bN1NWVsaIESN49tlniY+Pp23btoEO0ZigZQnCNEo7d+50TXtx+vRpunbtyuzZs7npppsqdM4bY/ynWSQIVbURLH7iyz6skydP8t5775GYmMju3buJiIhg4sSJ3HTTTYwZM8bqNxvTwII+QbRs2ZITJ07QuXNnSxI+pqqcOHHCdWNefZSUlLB582YSExP56KOPKC4uZsiQITz11FNMnTqVjh07+i5gY0ydBH2CiIqKIjMzk+PHjwc6lKDUsmVLoqKi6rzdvn37ePfdd1m7di3Hjx+nU6dO3HHHHdx0000BG1psjKko6BNEeHi4a0oLE1i5ubm8//77JCYm8t133xEWFsY111zDjBkzuOaaawgPDw90iMYYN0GfIExglZWV8fnnn/Puu++yceNGzp8/z4ABA3jsscf4yU9+YvWZjWnELEEYvzh06BCrV69mzZo1ZGVl0b59e26++WZmzJjBkCFDrD/ImCbAEoTxmXPnzpGSkkJiYiJfffUVIsKYMWP4zW9+w4QJE1yTCxpjmoYaE4SIjASuBnoABcAO4ENVPenn2EwjM2bMGHJycqos79ChAxMmTCA5OZn8/Hz69OnDI488wo033ki3bt0CEKkxxhe8JggRuQv4d+AAsA1IA1oCY4HficgO4PeqeqgB4jSNgKfkAHDmzBmSk5O54YYbmDFjBiNGjLBLSMYEgepaEG2AMapa4OmXIjIMiAUsQRhSU1Np06ZNoMMwxviQ1wShqn+tbkNV3e7zaEyjpKo11qO25GBM8Kn19Jci8hMR+UJEtovIA7VYv5eIfCQiu0TkBxF52Lm8k4hsFJF05+NFXraPF5E0EdkrIo/V/p9kfEVV2bx5M7fffju/+MUvAh2OMaaBeU0QIlK5mPGdwJXA5cDcWuy7BJivqgOd2z0oIoOAx4BNqhoLbHK+rvzeocBfgRuAQcBtzm1NA1BVNm3axC233MIvf/lLsrOzeeqppwIdljGmgVXXB/GAOHoan1LVI8Bh4FmgDMiqaceqmg1kO5+fFZFdQE8gARjvXO0N4GPgd5U2Hw3sVdX9ACLytnO7nbX6V5l6KSsrIyUlhRdffJG0tDSioqL4wx/+wI033khERAR/+9vfPHZU281uxgSn6vog7ne2Il4Wka3A74F/A1oDf6jLm4hINDAc+ALo6kweqGq2iHiqIN8TR0Iqlwlc4WXfs4HZAL17965LWMappKSE999/n5dffpl9+/bRt29fnnvuOaZNm0ZY2L8+Ilu2bAlglMaYhlbtfRCq+i2QICI/AdYCb6jq3+vyBiLSFngX+JWq5tZy+KOnlTzOK62qS4Gl4Cg5WpfYmruioiKSkpJYunQphw4dYsCAAbzwwgtMnjzZptY2xlR7H8Qc4H4cJ+bFQDyOy04pwP+nqptr2rmIhONIDm+paqJz8VER6e5sPXQHjnnYNBPo5fY6ilpc1jK1c/78ed59911eeeUVsrKyGDx4MH/961+57rrrrGynMcal2j4IVR0qIhHAZ6r6NvCfIvJ3HJebqk0Qzv6LV4Fdqvpnt1+tBX4BLHI+JnnY/CsgVkT6Aj8Cs4Dba/lvMl4UFBTwzjvv8Oqrr3Ls2DGGDx/O008/zbhx4+zGNmNMFdUliB9F5A9AK2B3+UJVPQU8Wot9j8Ex8ul7EdnuXPY4jsSwSkTuxXGT3S0AItIDWKaqU1S1RETmASlAKPCaqv5Qp3+ZccnLy2PFihX893//NydPnmT06NEsXryYK6+80hKDMcYr8VYy0tlymAwUAxtVtbQhA6uPkSNH6tatWwMdRqORm5vL8uXLWb58OWfOnGHs2LHMnTuXkSNHBjo0Y0wjISLbVNXjSaG6FkQPVX2vmp0K0FNVMy80QONbJ0+e5I033uDNN98kLy+P6667jrlz5zJ06NBAh2aMaUKqSxB/EpEQHH0E24DjOCbr6w9cC0wAFuLoUDaNwPHjx3nttddYuXIlhYWFTJo0iblz5zJw4MBAh2aMaYKquw/iFufdyz8D7gG6A/nALmA98KyqFjZIlKZa2dnZvPrqq6xatYri4mKmTp3KnDlz6N+/f6BDM8Y0YTXdB7ETeKKBYjF1dPjwYV555RUSExNRVRISErj//vvp06dPoEMzxgQBqyjXBB04cICXX36ZtWvXEhISwk9/+lPuu+8+evbsGejQjDFBxBJEE7Jnzx5eeuklPvjgAyIiIrjjjju499576dq1a6BDM8YEIUsQTcAPP/zAiy++yMaNG2ndujX33nsvd999N507dw50aMaYIFabmtTvAq8BH6hqmf9DMuW2b9/Oiy++yMcff0y7du148MEHufPOO7noIo8lNIwxxqdq04J4EbgbxzQb/wO8rqq7a9jGXIAvv/ySF198kU8//ZSOHTvyyCOP8LOf/Yx27doFOjRjTDNSY4JQ1Q+BD0WkA3AbsFFEDgOvAG+qarGfY2wWVJVPP/2Uv/3tb2zdupXIyEh++9vfMmvWLCvnaYwJiFr1QYhIZ+AOHHMrfQO8BYzFMdneeH8F1xyoKh9//DEvvvgi3377LV27duXJJ5/klltuoWXLloEOz5gGl5aWxuuvv85dd91FXFxcoMNp1mqc21lEEnHM3Noa+ImqTlfVd1T1IaCtvwMMVuXV22666SbmzJlDTk4OzzzzDB9++CF33nmnJYcgkpaWxoIFC0hLSwt0KI1eWloaS5Ys4eTJkyxZssSOWYBVN1lfec2G61T1Hw0cV700hcn6SktLWb9+PS+99BJ79+4lOjqaOXPmMG3aNMLDwwMdnvGx8hNeUVERERERzJs3z74Ve+F+rMrZMfO/6ibrqy5BfABchKNmdDKQqqol/grSFxpzgiguLmbt2rUsXbqUjIwMYmNjmTt3LvHx8Va9LUjV94SnqpSVlVV59LSspseatr/Qffpq36dPn2b37t2UlVUdKBkSEsLIkSPp1q0bYWFhhIeHExER4Xru7ScsLIyIiAjX82AthnWhl+TqlSCcG7bE0cdwA476DodwJItkVT1U50j8rDEmiKKiIhITE1m6dCk//vgjgwYNYu7cuUycODFoP7DBSlU5f/48+fn5FBQUVPt47Ngx9u7di7e/rxYtWiAiHk+Y1f1NNiUiQkhIiOvR/XnlZWfOnPGYHHwpNDS0xoRS22RT1yTlr791X7RQ650gPOyoL45kEQ90U9XRdYrEzxpTgigoKOB//ud/WLZsGUePHuWyyy7jgQce4JprrrEiPQFSVlbmOsGXn8wrn9grn+wr/76mv5eIiAhat25Nbm5utSe8Fi1aMHbs2FqdRGt6rO1JuL77rO37uD+W/9SWp9aW+zGdN28e/fv3p7i4mJKSEoqLi73+lJSUUFRUVGHdoqKiareraZ+lpRdWDickJKROyac26x07doz169dTUvKvCzv1SRI+SRAi0p6Ko57yVLXq/2YANYYEkZeXx9tvv81rr73GiRMnGDVqFA888ABXXXVV0CSGQI0yKSsrq3LS9nZi93SCLygoqPEE36JFC1q3bk2rVq1cj+7PW7duXe3vw8IcfyK1OeHZdfWKGnMfRFlZWa2TjbfkU5sk5Wkd9wRQG3U9ZheUIETkfuAZoAAoX1lVtV/tQ24YgUwQubm5vPnmm7zxxhucPn2aMWPGMHfuXEaNGhWQePzlQpq0paWlrhN1XU7s5esVFtY8u3zLli1dJ3H355VP8p4eW7Vq5dP+oMZ8wmusrFO/KvfkVP64ePFizpw543WbTp068cc//rFW+7/QBJEOXKWqObV6twAKRII4deoUb7zxBn//+9/Jy8vj2muvZe7cuVx22WUNGkdDSEtL47/+678oLv7XvZFhYWFce+21tG/f3uvJvfzx/Pnz1e5fRFwn6ppO6J6+vbdq1arR9evYCa/u7D6ImvmyhXqhCSIZmKGq+bV6twBqyASRk5Pjqt6Wn5/P5MmTmTNnDoMGDWqQ9/clVSU/P5/c3FzOnDnj8ef48eOcOnWq2v2IiNfLMtWd2MuXt2jRotGd4H3BTnjGH3zVQr3QBDEc+G/gC8D1FVBV/72G7V4DpgHHVHWIc9k7QHnkHYHTqjrMw7YZwFmgFCjxFnxlDZEgjh49yiuvvOKq3jZlyhTmzJlDbGysX9+3PsrKyjh79myVk33lRJCbm1uhVVAuPDycDh060KFDBw4fPuzx20q5iy66iD/+8Y9B089iTFPg71FMtZlq42XgH8D3QF3Gob0OLAGWly9Q1VvdgvoPwPtFNLi2MV3WyszMZOnSpa7qbdOnT2f27Nn07du3wWMpLi6u9oRf/nP27FmPnbKtW7emQ4cOtG/fnpiYGFcSqPzTsmVL1wm/pibt3XffbcnBmAYWFxfHvHnz/NZCrU0L4lNV/bd67VwkGlhX3oJwWy447qm4TlXTPWyXAYysa4LwRwsiIyPDVb1NRLj55pu57777iIqK8un7qCqFhYVeL/G4J4L8/KpX+0SEdu3aeTzRt2/fvsLr+t6xbZ2uxgSfC21BfCQis4H3qHiJ6eQFxHQ1cNRTcijfPbBBRBR4WVWXetuRM7bZAL17965XMJ6uEaenp/PSSy+xfv16wsPDuf3227n33nvp1q1bnfZdVlZGXl5ejZd5zpw54/EyT1hYmOvE3q1bN+Li4qqc8Dt06EC7du38fv2+/NuKdboa0zzUpgVxwMPiWg1zraYF8SKwV1X/w8t2PVQ1S0QuBjYCD6nqJzW9X11bEGPGjKGkpIRLL72U0NBQSktL+f7778nPz6eoqIjWrVtz2223cc899xAZGVlh25KSkmo7dcuTgLcbplq1auXx233lb/6tW7dudJdurNPVmODhszup6/HG0VRKECISBvwIjFDVzFrs42kcN+U9X9O6dU0QV1xxhSs5lCstLWXPnj3Ex8fzb//2b5SUlHj8xn/u3DlPsdK2bdtaXeaJiIiodZzGGOMvF3SJyTkf0wM46j8ojqm/X1LVmu9a8mwisNtbchCRNkCIqp51Pp+E40Y9n0pLS6uSHMAxX8vAgQM5ePAgBw8eBByXecpP8BdffDGxsbFeL/PYxHvGmGBRmz6I5TiGnP6X8/VtwN+BW6rbSERW4pjoL1JEMoGFqvoqMAtYWWndHsAyVZ0CdAVWOy+rhAErVDW5tv+g2nr99derPZm3a9eORx99lPbt29OmTZtGd5nHGGP8rTZ9EN+q6mU1LWsM6nKJKS0tjT/96U8ek0RpaSm/+c1v7Pq6MSboVXeJqTbDXr4RkSvddnYFsMVXwQVKXFwc33//fZVZGss7qi05GGOau9okiCuAT0Ukw3l/wmfANSLyvYh859fo/CwsLKxCkihPDuUzchpjTHNWmzNhvN+jCJAtWxwNIRu2aYwxVfl1mGtDawz1IIwxpim50DupjTHGNEKpqamsWrWKnJwcIiMjmTlzJmPHjvXZ/i1BGGNME5SamsqyZctcc6Pl5OSwbNkyAJ8lieCbfN8YY5qBVatWVZlduaioiFWrVvnsPawFYYwxTYCqkpOTQ1paGmlpaeTkeJ7s+sSJEz57T0sQxhjTCJWVlXHo0CFXQtizZw8nTzom0W7VqhXh4eEeZ4Du3Lmzz2KwBGGMMY1AYWEh+/btcyWEvXv3UlBQAECnTp2Ii4tz/fTq1YtPP/20Qh8EOOqzzJw502cxWYIwxpgAOH36NHv27HElhIyMDMrKyhARoqKiGDNmjCshVC43AP/qiF61ahUnTpygc+fOPh/FZPdBGGOMn6kqWVlZFRLC0aNHAUft95iYGFcyiI2NpU2bNg0Wm90HYYwxDaikpIT9+/e7+g727NnD2bNnAcdM0QMGDGDChAnExcXRt2/fRju9T+OMyhhjmpBz5865Wgd79uxh3759rg7kbt26MXz4cFcLoXv37k2mfIAlCGOMqYPKw0337NlDZmYmqkpoaCjR0dFMnDjRlRA6dOgQ6JDrzRKEMcZUo6bhprGxsVx55ZUMGDCAmJgYWrZsGeCIfccShDHGuKnrcNOQkOCdkMIShDGmWTtz5owrGaSlpXHw4EFKS0trPdw0mFmCMMY0G7UZbjpt2rSADDdtjPyWIETkNWAacExVhziXPQ3cBxx3rva4qq73sG088BcgFFimqov8FacxJniVDzd1H2HUFIebBoo/j8brwBJgeaXlL6jq8942EpFQ4K/A9UAm8JWIrFXVnf4K1BgTHMqHm5YnhGAZbhoofksQqvqJiETXY9PRwF5V3Q8gIm8DCYAlCGOMi/tw0/KEEKzDTQMlEO2peSLyc2ArMF9VT1X6fU/gsNvrTOCKhgrOGBM41VVIa87DTQOloRPEi8AfAHU+/gdwT6V1PLX5vE4YJSKzgdkAvXv39k2UxpgG56lC2tKlS9m6dSsFBQXNerhpoDRoglDVo+XPReQVYJ2H1TKBXm6vo4Csava5FFgKjsn6fBOpMaYhqSorV66sUiGtpKSEL7/8kl69ejXr4aaB0qAJQkS6q2q28+VNwA4Pq30FxIpIX+BHYBZwewOFaIxpAEVFRWRkZJCenu76OXWq8tXmf3nuuecaMDpTzp/DXFcC44FIEckEFgLjRWQYjktGGcD9znV74BjOOkVVS0RkHpCCY5jra6r6g7/iNMb434kTJyokg4yMDEpKSgDo0qULAwcO5Ntvv+XcuXNVtrXWQuD4cxTTbR4Wv+pl3Sxgitvr9UCV+yOMMY1f5dbB3r17XZ3JERER9OvXjxtuuIH+/fsTGxtLx44dgap9EOXr+7JCmqkbuyvEGHNBamodXHLJJcTGxhIbG0vv3r293ozWEBXSTN1YRTljTK1V1zoon6oiNja2SuvANF5WUc4YUy++ah2Ypsn+N40xABQXF3PgwIFqWwee+g5M8LIEYUwzZa0DUxP7HzemGaipddCvXz/i4+NdCcFaBwYsQRgTlKx1YHzBPhXGNHHWOjD+YgnCmCamNq2D8o7kPn36WOvA1Jt9coxpxNxbB3v37iU9Pd1r66B///5cdNFFAY7YBBNLEMY0IjW1DsprJVvrwDQE+3QZ40fVFcApLi6uMqOptQ5MY2IJwhg/8VYAJzU1lfz8/Aqtg8jISGsdmEbHPoHG+EFBQQFvvfWWxwI43333HXFxccTHx7s6k611YBojSxDGXKCSkhIOHTrEvn37XD9ZWVlUNxHmwoULGzBCY+rHEoQxdVBWVkZ2drYrEezfv5+DBw+6LhW1b9+emJgYrrrqKjZs2EBubm6VfVgBHNNUWIIwxgtV5eTJkxVaBgcOHKCgoACAli1bujqSY2Ji6NevH5GRkYgIABdffLEVwDFNmiUIY5zy8vJcrYLyx9OnTwMQGhpKnz59GDt2LP369SMmJoYePXoQEhLidX9WAMc0dZYgTLN0/vx5MjIyKiSEo0ePAiAi9OjRg0svvdTVMujTpw/h4eF1fp+xY8daQjBNliUIE/RKS0s5fPhwhZbB4cOHKSsrA6Bz587ExMRw7bXXEhMTQ9++fWndunWAozYm8PyWIETkNWAacExVhziX/Qn4CVAE7APuVtXTHrbNAM4CpUCJt3J4xlSmqhw9erRCyyAjI8PVD9CmTRtiYmIYPnw4MTExxMTE2OR1xnjhzxbE68ASYLnbso3AAlUtEZHngAXA77xsf62q5vgxPhMETp065UoE5Unh3LlzgKNDODo6mgkTJriSwcUXX+zqRDbGVM9vCUJVPxGR6ErLNri9/Bz4qb/e3wSf/Px8Dhw4UGFUUfnUFCEhIfTq1YvRo0fTv39/+vXrR1RUFKGhoQGO2pimK5B9EPcA73j5nQIbRESBl1V1qbediMhsYDZA7969fR6kCYyioqIKN5/t37+frKws1++7du3KJZdc4moZ9OnThxYtWgQwYmOCT0AShIg8AZQAb3lZZYyqZonIxcBGEdmtqp94WtGZPJYCjBw50vutq6bRKisrIysri71797ouFx06dIjS0lIAOnbsSExMjGuIab9+/Wjbtm2AozYm+DV4ghCRX+DovJ6gXuYiUNUs5+MxEVkNjAY8JgjTtKgqOTk5FVoGBw4coLCwEIBWrVrRr18/pk6d6rrfoFOnTtZvYEwANGiCEJF4HJ3S16hqvpd12gAhqnrW+XwS8EwDhml8KDc3t8rNZ+XTT4SFhREdHc24ceNc9xt079692pvPjDENx5/DXFcC44FIEckEFuIYtdQCx2UjgM9VdY6I9ACWqeoUoCuw2vn7MGCFqib7K05TN9XVNygsLOTAgQMVRhUdP34ccNx81rNnT4YPH+5qGfTu3dumtDamEZPqZpxsakaOHKlbt24NdBhBq3J9A3BMQTFgwADy8vLIzMx0zWDapUsXVyKIiYkhOjqaVq1aBSp0Y4wXIrLN271m9vXN1EphYSFvvvlmlfoGpaWl7N69m8suu4xRo0a5LhV16NAhQJEaY3zFEoTxqLCwkPT0dHbu3MmuXbvYt2+fa1RRZarKb3/72waO0Bjjb5YgDOA9IYSGhrpGFX388cdW38CYZsQSRDNVWFjInj17XAlh//79VRLCoEGDGDBgAC1btgQgKirK6hsY04xYgmgm6pMQKrP6BsY0LzaKKUjVlBAGDhxYY0IwxgQ/G8XUDFSXEGJiYpg2bRqDBg0iNjbWEoIxplYsQTRRlhCMMf5mCaKJKCwsJC0tjV27drFz504OHDhgCcEY41eWIBopSwjGmECzBNFIVE4I+/fvp6yszBJCEDh//jw7d+5k0KBBVrOiFux4NR6WIAKkpoQwffp0Bg4c2CgTgv0B101GRgZnzpwhIyODuLi4QIfT6NnxajwsQTSQgoKCKp3KTSUhVNbY/oBVFVWlrKzMJ899ta/Kd51nZ2eTnZ0NQOvWrT3WuKi8rCmuU9/tyo+N++vs7GxEhOjoaESEkJCQah9rs477usFQZ8SfX9gsQfhJMCUEcFR927x5M+73zbj/AcfFxQXsxNwQ3E8onk5GlZ+HhobSsWNH8vPzK9x53qJFC9q0aVOh5oWne5Fqs6y6dSo/1nc/DblOWFgYpaWlHrc5cOBAle184UKTjC8SVV229cSfX9gsQfhIU0kIpaWllJSUUFxcXOXR07Lyx+pOxKrK7t27q33f2p5Y3f8gwsLC6rRNfZ7XZb36SEtLIzs7m5CQEMrKyujUqVOjaHU1VpWPV/fu3RkwYIDXLwn1fazPNuXJq6Z9+JP730dxcXGF35V/YQsJCWHcuHE+eT9LEPVUXULo378/06dPd3Uq++M6vT9O9CJCeHg4YWFhhIeH07JlS8LDwyssO378OKdOnUJEUFW6dOlCdHR0jSfZ5qq4uJgePXrQo0cPsrKyqkyXbirydLzcP0tNgb8TVfljcXExubm5nD9/HnAkj8jISGJiYnz2b7EEUUv+SgiBOtFXfgwPD6/VyfzkyZNV/oDbtGlT639vczNkyBDX8wEDBgQwkqYhGI5X+SXGhlC5xRUaGurTL6TNPkGkpqaydu1ahg0bxjfffENCQgJjx46loKCgyn0I1SUE9xN9QUGBX0/07id1Tyf80NBQv31rD4Y/YGOChb9bqM16sr7yEpqjRo1i0KBB7Ny5k6+++opu3bpx7tw5IiIiaNWqFb169aJnz5506dKFdu3aUVZWdkEnek/f3hv6RG+MMWCT9Xl1/vx57rrrLtfrwYMHM3jwYK/rFxYWUlpa2qi+0RtjjL/4LUGIyGvANOCYqg5xLusEvANEAxnATFU95WHbeOAvQCiwTFUX+SPGFStWcOWVVxITE+O6hnfy5En27t3LPffcYyd6Y0yz5s9hAa8D8ZWWPQZsUtVYYJPzdQUiEgr8FbgBGATcJiKD/BFgmzZtKC4uRkQoKSlBRDh69ChZWVlcfPHFXHTRRbRr146WLVsSFhZmycEY06z4LUGo6ifAyUqLE4A3nM/fAG70sOloYK+q7lfVIuBt53Y+N3PmTNq0acPOnTtZs2YNO3fupG3btlZC0xhj8G8LwpOuqpoN4Hy82MM6PYHDbq8zncs8EpHZIrJVRLYeP368TsGMHTuWmJgY0tLSOHXqFGlpacTExFgJTWOMoXF2Unu6juN1qJWqLgWWgmMUU13fbOzYsZYQjDHGg4ZuQRwVke4AzsdjHtbJBHq5vY4CshogNmOMMW4aOkGsBX7hfP4LIMnDOl8BsSLSV0QigFnO7YwxxjQgvyUIEVkJfAbEiUimiNwLLAKuF5F04Hrna0Skh4isB1DVEmAekALsAlap6g/+itMYY4xnfuuDUNXbvPxqgod1s4Apbq/XA+v9FJoxxphaaBrTIxpjjGlwQTUXk4gcBw7Wc/NIIMeH4QQ7O151Y8erbux41c2FHK8+qtrF0y+CKkFcCBHZ6m3CKlOVHa+6seNVN3a86sZfx8suMRljjPHIEoQxxhiPLEH8y9JAB9DE2PGqGztedWPHq278crysD8IYY4xH1oIwxhjjkSUIY4wxHgVtghCRXiLykYjsEpEfRORh5/JOIrJRRNKdjxc5l3d2rp8nIkvc9tNORLa7/eSIyP8N0D/Lb+pxvEa7HZNvReQmt32NEJHvRWSviPynBGGlpboeL7ftejs/Y792W2bHq+rnK1pECtw+Yy+57cuOl4fPl4gMFZHPnOt/LyItncvrf7xUNSh/gO7A5c7n7YA9OCrULQYecy5/DHjO+bwNMBaYAyypZr/bgHGB/vc1guPVGghz2/aY2+svgatwTN3+AXBDoP99gT5ebtu9C/wP8Gu3ZXa8qn6+ooEdXvZlx6vq8QoDvgMuc77uDIRe6PEK2haEqmar6tfO52dxTPzXEy9V7VT1nKqmAoXe9ikisTiKHG32X+SBUY/jla+OiRUBWuKs2eGcxr29qn6mjk/ncjxXDmzS6nq8AETkRmA/8IPbMjteNVeZdLHj5fV4TQK+U9VvnducUNXSCz1eQZsg3IlINDAc+ILaVbXz5jbgHeeBDlq1PV4icoWI/AB8D8xxJoyeOGp6lKu2ImAwqM3xEpE2wO+A/1Npczte3v8e+4rINyLyTxG52rnMjpfn4zUAUBFJEZGvReS3zuUXdLwaY0U5nxKRtjia9b9S1dwLvFw5C7jTJ4E1UnU5Xqr6BTBYRAYCb4jIB9SxImBTV4fj9X+AF1Q1r9I6drw8ywZ6q+oJERkBrBGRwdjx8rZqGI5L5KOAfGCTiGwDcj2sW+vjFdQJQkTCcRzct1Q10bn4qIh0V9Vs8V7VztO+LsNxjX2bn8INuPoeL1XdJSLngCE4vqFEuf06aCsC1vF4XQH8VEQWAx2BMhEpdG5vx6vS8VLV88B55/NtIrIPx7dk+3x5/nxlAv9U1RzntuuBy4E3uYDjFbSXmJw99a8Cu1T1z26/qk1VO09uA1b6LsLGpa7HSxwV/8Kcz/sAcUCGs9l7VkSudO7z59T+GDcZdT1eqnq1qkarajTwf4H/X1WX2PHy+vnqIiKhzuf9gFhgvx0vr+evFGCoiLR2/l1eA+y84OPV0L3zDfWDo7mlOHr2tzt/puDo3d8EpDsfO7ltkwGcBPJwZORBbr/bD1wS6H9XYzleOC61/eBc72vgRrd9jQR2APuAJTjv2A+mn/p8vty2fZqKo5jseFX9fN3s/Hx96/x8/cSOV43nrzucx2wHsNgXx8um2jDGGONR0F5iMsYYc2EsQRhjjPHIEoQxxhiPLEEYY4zxyBKEMcYYjyxBGGOM8cgShDHGGI/+H68L2kXNLbseAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "pmoysf       = np.zeros((2,4))\n",
    "pmoys0f      = np.zeros((2,4))\n",
    "pmoysf[0,:]  = yy_eu\n",
    "pmoysf[1,:]  = yy_us\n",
    "pmoys0f[0,:] = yy_eu0\n",
    "pmoys0f[1,:] = yy_us0\n",
    "\n",
    "\n",
    "xx      = np.zeros((2,5))\n",
    "xx[0,0] = PMoY10[0]\n",
    "xx[1,0] = PMoY10[-1]\n",
    "zz      = np.zeros((2,5))\n",
    "zz[0,0] = PMoY10[0]\n",
    "zz[1,0] = PMoY10[-1]\n",
    "for jj in range(2):\n",
    "    for ii in range(horizon):\n",
    "        xx[jj,ii+1] = pmoysf[jj,ii] \n",
    "        zz[jj,ii+1] = pmoys0f[jj,ii]\n",
    "    \n",
    "lyhorizon = ('2017','2030','2040','2050','2060')\n",
    "\n",
    "plt.plot(np.arange(len(lyhorizon)),list(xx[0,:]),'o-',label=r'EU, $\\Delta p \\neq 0$',color = '0.35')\n",
    "plt.plot(np.arange(len(lyhorizon)),list(xx[1,:]),'s-',label=r'US, $\\Delta p \\neq 0$',color = '0.15')\n",
    "plt.plot(np.arange(len(lyhorizon)),list(zz[0,:]),'*-',label=r'EU, $\\Delta p = 0$',color = '0.75')\n",
    "plt.plot(np.arange(len(lyhorizon)),list(zz[1,:]),'D-',label=r'US, $\\Delta p = 0$',color = '0.4')\n",
    "plt.xticks(np.arange(len(lyhorizon)), labels=lyhorizon)\n",
    "plt.legend() #(loc='center left', bbox_to_anchor=(1, 0.5))\n",
    "plt.ylabel('pm/y (%)')\n",
    "plt.savefig('forecasts.eps',dpi=600)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b73ef8d1",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
